###########################################################################
###########################################################################
####Non-state Atrocities in Capital Cities - Mechanism 1 (Regressions) ####
###########################################################################
###########################################################################
library(MASS) 
library(pscl)
library(foreign)
library(stargazer)

# Set working library
setwd("~/Data/Mechanism 1/")

#################
###Afghanistan###
#################
# Read in data for analysis
afg.analysis <- read.dta("afganal.dta")

###Baseline model
afg.dat.1 <- na.omit(afg.analysis[,c("incidentnonstatefull", "cit_ind", "capital", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03", "year04", "year05", "year06", "year07", "year08")])
#Run model
afg.nb.1 <- glm.nb(incidentnonstatefull ~ capital + cit_ind + year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03 + year04 + year05 + year06 + year07 + year08,
                   data=afg.dat.1)
summary(afg.nb.1)
AIC(afg.nb.1)

###Full model
afg.dat.2 <- na.omit(afg.analysis[,c("incidentnonstatefull", "cit_ind", "lagincidentnonstatefull", "loglagpop",  "lagcivconflagtemp", "capital", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03", "year04", "year05", "year06", "year07", "year08")])
#Run model
afg.nb.2 <- glm.nb(incidentnonstatefull ~ capital + cit_ind + lagincidentnonstatefull +
                     lagcivconflagtemp + loglagpop + 
                     year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03 + year04 + year05 + year06 + year07 + year08,
                   data=afg.dat.2)
summary(afg.nb.2)
AIC(afg.nb.2)

##############
###Pakistan###
##############
# Read in data for analysis
pak.analysis <- read.dta("pakanal.dta")

###Baseline model
pak.dat.1 <- na.omit(pak.analysis[,c("incidentnonstatefull", "cit_ind", "capital", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03", "year04", "year05", "year06", "year07", "year08")])
#Run model
pak.nb.1 <- glm.nb(incidentnonstatefull ~ capital + cit_ind + year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03 + year04 + year05 + year06 + year07 + year08,
                   data=pak.dat.1)
summary(pak.nb.1)
AIC(pak.nb.1)

###Full model
pak.dat.2 <- na.omit(pak.analysis[,c("incidentnonstatefull", "cit_ind", "lagincidentnonstatefull", "loglagpop",  "lagcivconflagtemp", "capital", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03", "year04", "year05", "year06", "year07", "year08")])
#Run model
pak.nb.2 <- glm.nb(incidentnonstatefull ~ capital + cit_ind + lagincidentnonstatefull +
                     lagcivconflagtemp + loglagpop + 
                     year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03 + year04 + year05 + year06 + year07 + year08,
                   data=pak.dat.2)
summary(pak.nb.2)
AIC(pak.nb.2)

##########
###Iraq###
##########
# Read in data for analysis
irq.analysis <- read.dta("irqanal.dta")

###Baseline model
irq.dat.1 <- na.omit(irq.analysis[,c("incidentnonstatefull", "cit_ind", "capital", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03", "year04", "year05", "year06", "year07", "year08")])
#Run model
irq.nb.1 <- glm.nb(incidentnonstatefull ~ capital + cit_ind + year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03 + year04 + year05 + year06 + year07 + year08,
                   data=irq.dat.1)
summary(irq.nb.1)
AIC(irq.nb.1)

###Full model
irq.dat.2 <- na.omit(irq.analysis[,c("incidentnonstatefull", "cit_ind", "lagincidentnonstatefull", "loglagpop",  "lagcivconflagtemp", "capital", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03", "year04", "year05", "year06", "year07", "year08")])
#Run model
irq.nb.2 <- glm.nb(incidentnonstatefull ~ capital + cit_ind + lagincidentnonstatefull +
                     lagcivconflagtemp + loglagpop + 
                     year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03 + year04 + year05 + year06 + year07 + year08,
                   data=irq.dat.2)
summary(irq.nb.2)
AIC(irq.nb.2)


############
###India###
###########
ind.analysis <- read.dta("indanal.dta")

###Baseline model
ind.dat.1 <- na.omit(ind.analysis[,c("incidentnonstatefull", "cit_ind", "capital", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03", "year04", "year05", "year06", "year07", "year08")])
#Run model
ind.nb.1 <- glm.nb(incidentnonstatefull ~ capital + cit_ind + year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03 + year04 + year05 + year06 + year07 + year08,
                   data=ind.dat.1)
summary(ind.nb.1)
AIC(ind.nb.1)

###Full model
ind.dat.2 <- na.omit(ind.analysis[,c("incidentnonstatefull", "cit_ind", "lagincidentnonstatefull", "loglagpop",  "lagcivconflagtemp", "capital", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03", "year04", "year05", "year06", "year07", "year08")])
#Run model
ind.nb.2 <- glm.nb(incidentnonstatefull ~ capital + cit_ind + lagincidentnonstatefull +
                     lagcivconflagtemp + loglagpop + 
                     year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03 + year04 + year05 + year06 + year07 + year08,
                   data=ind.dat.2)
summary(ind.nb.2)
AIC(ind.nb.2)


######Export to LaTex
stargazer(afg.nb.1, afg.nb.2, pak.nb.1, pak.nb.2, ind.nb.1, ind.nb.2, irq.nb.1, irq.nb.2)